Large morphological datasets

Plantkon imagers can image planktonic organisms from 300nm (FlowCam nano) up to 10cm (ISIIS) (Lombard et al, 2019). Some run in a lab but more and more are automated and can be deployed on cruises or on fixed platforms, generating millions and millions of images. Those datasets are mostly used for taxonomic puposes. However, a number of morpgological characteristics associated to each individual are often generated to be used within machine learning algorithms. If those features form the basis of the plankton image recognition, they are also interesting measures of variability within the community that do not depend on the species level.

Such morphological datasets that combine numerous columns of morphological features and millions of observations, can’t easily be opened in any of the widely used programming softwares (Python, R, MATLAB) for further exploration with multivariate data analysis. Even if it was possible, running a Principal Component Analysis or a t-SNE on millions of lines takes time and computer power. Knowledge of methods and ressources specific to big data analysis such as MapReduce with Hadoop or High Performance Computing could be an answer but the ressources are not available to everyone and the tools have a quite steep learning curve.

A methodology based on subsampling, dimension reduction and clustering

Here, we detail a workflow to explore morphological variation and diversity of planktonic organisms based on subsampling, Principal Component Analysis, K-means clustering and reprojection (Figure 1). The goal is to explore the main morphological variations within the dataset with a dimension reduction technique, the PCA, and to cluster the images projected in this morphological space into morpho groups that can be used similarly as species with the diversity indices developed by Villeger et al. (2008). We use subsampling to run these analysis and the single value decomposition (SVD) matrix from the PCA with the k nearest neighbour method to reproject the totality of the dataset into the morphological space and assign them a morpho group. The methodology comports 3 steps:

1. How many images have to be subsampled to capture the morphological variation?
2. How many clusters are needed to capture the morphological diversity?
3. Analysis and reprojection

Figure 1. Workflow of the methodology based on subsampling and a combination of dimension reduction, clustering and diversity indices

Applications and example

The full methodology was applied to a 2-years dataset of phytoplankton imagery from an Imaging FlowCytobot (Sonnet et al., in preparation) while the morphological approach combining PCA, k-means and diversity and designed by Jean-Olivier Irisson and Sakina-Dorothée Ayata was used on Mediterannean zooplankton (Irisson et al., pre-print) Arctic zooplankton (Vilgrain et al., 2020, in press) and marine detritus (Trudnowska et al., 2020, in review).

Below, you can find an example with a morphological dataset from an Imaging FlowCytobot. The graphs come from the analysis of the 2-years dataset but if you want to try running the analysis, in the github folder you will find a 2-days extract of the dataset and the files that stem from this analysis applied on it.

Storage

I assume the morphological features are stored for each image within a database. All the queries within this document are only examples and will need to be adapted to the structure of your database. The same applies for the database information and password. For reference, Figure 2 shows the structure of the database I used (designed by Audrey Ciochetto).

If you only have csv files, you can combine them within the shell and work from there and the principales of the analysis are the same. I also have some codes that I previously used when I was in this situation and I would be happy to share them if they could be useful.

Figure 2. Diagram of the database structure in the Mouw lab. Main boxes are tables and each table has a first column called id. Smaller boxes on arrows are the column name that relates the right table with the id column of the left table.

Codes

Most of the codes have not been turned into functions because there are many possibilities in your column names, morphological characteristics, conditions specific to the deployment of any plankton imager and the parameters you will want to test. However, both within this document and at the beginning of each code, there is an indication of which sections of the code you might be willing to adapt.

Throughout the document and the codes associated, I assume you are in an R project and directory that have 4 folders: data, figures, documents and scripts.

1.1. Random files

Select a percentage of random samples from your dataset: R script 1-1_sample_files.R.

1. Get the number of files present in the database and meeting your criteria (you might want to change these criteria in the query)
2. Calculate the number of files to sample based on a percentage (the percentage is 5 in the code)
3. Query the files id and take randomly the number determined in 2. (you might want to change these criteria in the query)
4. Query the images that meet your criteria from the files id (you might want to change the criteria of the images in the query)
5. Extract counts per sample
6. Export to csv

Outputs: combined images file (\(\small{\textit{subsp_combined.csv}}\)) and counts per file (\(\small{\textit{subsp_counts_cleaned.csv}}\))

Extract column names in bash:

head -1 data/subsp_combined.csv > data/colnames.csv

1.2. Subsampling

Subsample 3, 10, 100, 250, 500 and 1000 images per sample and repeat 10 times (total: 60 files): bash script 1-2_subsampling.sh

This code calls 2 other codes, extract.awk and indices.R.

Arguments:

• -d: path where to get and write the data
• -p: path where to get the scripts (no slash at the end)
• -f: filename of the combined images (default: combined.csv)
• -c: filename with the counts per sample after query/cleaning (default: counts_cleaned.csv)
• -s: type of sampling, proportion of the file (proportion) VS number of images per sample (date) (default)
• -i: replicate number
• -n: number of images to sample per file or proportion of the file to sample

Output: \(\small{\textit{combined_sampled{i}_{n}im.csv}}\)

# bash code
for i in 3 10 50 100 250 500 1000 # images subsampled by sample (you might want to try other ones)
do
    for j in {1..10} # replicates  
    do
      bash scripts/1-2_subsampling.sh -d data/ -p scripts -f subsp_combined.csv \
      -c subsp_counts_cleaned.csv -s date -i $j -n $i
    done
done
rm data/idxToSample.csv # remove the csv file with the random indices generated at each iteration  

Sometimes, a couple of the subsamples have a header. To make sure none of them have one, remove lines that contain the header pattern (make sure this pattern does not appear in your filename!).

grep "id" data/combined_sampled* | cut -d: -f1 | xargs -n1 sed -i .bak '1d' # select the files with a header and remove this header
rm data/*.bak # remove the temporary .bak files
grep "id" data/combined_sampled* # verify there are none left 

1.3. Subsample Analysis

mysql -h images.gso.uri.edu -D plankton_images -u student -pRestmonami -e "query;" > data/mincols.csv

SELECT MIN(Area) AS area, MIN(Biovolume) AS vol, MIN(ConvexArea) AS 'c.area', 
MIN(ConvexPerimeter) AS 'c.perim', MIN(Eccentricity) AS ecc, 
MIN(EquivDiameter) AS 'eq.diam', MIN(Extent) AS extent, MIN(H180) AS h180, MIN(H90) AS h90, 
MIN(Hflip) AS hflip, MIN(MajorAxisLength) AS 'maj.ax', 
MIN(MinorAxisLength) AS 'min.ax', MIN(Perimeter) AS perim, MIN(Solidity) AS solid, 
MIN(numBlobs) AS 'nb.blobs', MIN(shapehist_kurtosis_normEqD) AS 'kurt.dist', 
MIN(shapehist_mean_normEqD) AS 'mean.dist', MIN(shapehist_median_normEqD) AS 'med.dist', 
MIN(shapehist_mode_normEqD) AS 'mode.dist', MIN(shapehist_skewness_normEqD) AS 'skew.dist', 
MIN(summedArea) AS 's.area', MIN(summedBiovolume) AS 's.vol', MIN(summedConvexArea) AS 's.c.area', 
MIN(summedConvexPerimeter) AS 's.c.perim', MIN(summedMajorAxisLength) AS 's.maj.ax', 
MIN(summedMinorAxisLength) AS 's.min.ax', MIN(summedPerimeter) AS 's.perim', 
MIN(texture_average_contrast) AS 'tx.contrast', MIN(texture_average_gray_level) AS 'tx.gray', 
MIN(texture_entropy) AS 'tx.entropy', MIN(texture_smoothness) AS 'tx.smooth', 
MIN(texture_third_moment) AS 'tx.3rd.mom', MIN(texture_uniformity) AS 'tx.unif', 
MIN(RotatedBoundingBox_xwidth) AS 'r.bbox.x', MIN(RotatedBoundingBox_ywidth) AS 'r.bbox.y',
MIN(Area_over_PerimeterSquared) AS 'area.perim2', MIN(Area_over_Perimeter) AS 'area.perim', 
MIN(H90_over_Hflip) AS 'h90.hflip', MIN(H90_over_H180) AS 'h90.h180', 
MIN(summedConvexPerimeter_over_Perimeter) AS 's.c.perim.perim', 
MIN(MajorAxisLength/MinorAxisLength) AS 'maj.min', 
MIN(Perimeter/MajorAxisLength) AS 'perim.maj', 
MIN(4*3.14159265359*Area_over_PerimeterSquared) AS 'circ' FROM roi

find data/combined_sampled* ! -name combined_sampled1_250im.csv | xargs -n1 rm # change your sample name 

2.1. Compute clustering with different number of groups

On one of the samples for the number of images determined in 1, run: R script 2-1_explo_kmeans.R.

ATTENTION: this code calls the functions pcatrans.R, lib_functional_diversity.R and morphological_diversity.R. The function morphological_diversity is a wrap-up around the function lib_functional_diversity.R from Villéger that reearranges the dataset into weights per date per morph and centers coordinates necessary for the diversity indices.

1. Run the pca (with 4 components) on the subsample you have kept and store the coordinates for each image (you might want to change the number of components)
2. Perform the kmeans 10 times for 10, 25, 50, 100, 250, 500 and 1000 clusters => long step and you might want to change the numbers of clusters
3. Export the tibble with pca coordinates and cluster number, and the one with the kmeans centers
4. Compute the diversity for each of the output => longest step
5. Export the diversity per date

Outputs: the pca coordinates and cluster number for each image with different clusterings (\(\small{\textit{subsp_pca+conc+cluster.csv}}\)), the kmeans centers (\(\small{\textit{subsp_kmeans_centers.csv}}\)) and the diversity per date (\(\small{\textit{subsp_diversity_explo.csv}}\)).

2.2. Determine the number of groups to use

Use the csv files generated before in the R script 2-2_compare_kmeans.R. All parts are independant, if you only want to compute some of the graphs.

1. Make graphs based on the centers of the morpho groups: position of the centers for the different subsampling in the 1st and 2nd dimension (Figure 5) and in the 3rd and 4th dimension (Figure 6)

Figure 5. Position of the morpho groups’ centers in the first 2 axis of the morphological space created by the PCA, depending on the number of clusters chosen.

Figure 6. Position of the morpho groups’ centers in the third and fourth axis of the morphological space created by the PCA, depending on the number of clusters chosen.

2. Make graphs based on the number of morphs per days: correlation plot between the scenario of total clusters and time series (Figure 7)

Figure 7. (left) Correlation between the time series of morphological groups and (right) mean number of morphological groups per sample date for each number of clusters. The red curve is the smoothed mean number of morphological groups with a window of 10 values before and after and the error bars are in grey. See Appendix for all clusters.

3. Time series of diversity indices based on the total number of clusters (Figure 8)

Figure 8. Time series of diversity indices for 10, 25, 50, 150 and 1000 clusters.

4. Correlation plots of the time series of indices based on the total number of clusters (Figure 9)

Figure 9. (top) Correlation between diversity time series depending on the number of clusters and (bottom) Correlation between adjacent numbers of clusters: the rst dot is the correlation between the curves of 10 and 25 clusters, the second between 25 and 50 clusters, etc.

In my case, I have decided to settle for 150 clusters but as often in clustering, there is a piece of subjectivity involved.

You can delete the subsampled file.

# bash 
rm data/combined_sampled* 

3.1. Data

Same as for the subsampling, query the data through bash but this time query your whole dataset. .

mysql -h images.gso.uri.edu -D plankton_images -u student -pRestmonami -e  "query;" > data/combined.csv

For reference, here is the query I used and it took 411 minutes (7h) to pull 110 741 191 rows from the database (but I have very specific conditions in my query that make it much slower):

SELECT roi.id, filename, date, Area AS area, Biovolume AS vol, ConvexArea AS 'c.area', ConvexPerimeter AS 'c.perim', Eccentricity AS ecc, EquivDiameter AS 'eq.diam', Extent AS extent, H180 AS h180, H90 AS h90, Hflip AS hflip, MajorAxisLength AS 'maj.ax', MinorAxisLength AS 'min.ax', Perimeter AS perim, Solidity AS solid, numBlobs AS 'nb.blobs', shapehist_kurtosis_normEqD AS 'kurt.dist', shapehist_mean_normEqD AS 'mean.dist', shapehist_median_normEqD AS 'med.dist', shapehist_mode_normEqD AS 'mode.dist', shapehist_skewness_normEqD AS 'skew.dist', summedArea AS 's.area', summedBiovolume AS 's.vol', summedConvexArea AS 's.c.area', summedConvexPerimeter AS 's.c.perim', summedMajorAxisLength AS 's.maj.ax', summedMinorAxisLength AS 's.min.ax', summedPerimeter AS 's.perim', texture_average_contrast AS 'tx.contrast', texture_average_gray_level AS 'tx.gray', texture_entropy AS 'tx.entropy', texture_smoothness AS 'tx.smooth', texture_third_moment AS 'tx.3rd.mom', texture_uniformity AS 'tx.unif', RotatedBoundingBox_xwidth AS 'r.bbox.x', RotatedBoundingBox_ywidth AS 'r.bbox.y', Area_over_PerimeterSquared AS 'area.perim2', Area_over_Perimeter AS 'area.perim', H90_over_Hflip AS 'h90.hflip', H90_over_H180 AS 'h90.h180', summedConvexPerimeter_over_Perimeter AS 's.c.perim.perim', MajorAxisLength/MinorAxisLength AS 'maj.min', Perimeter/MajorAxisLength AS 'perim.maj', 4*3.14159265359*Area_over_PerimeterSquared AS 'circ', class_id 
FROM roi JOIN raw_files ON roi.raw_file_id = raw_files.id 
         JOIN auto_class ON roi.id = auto_class.roi_id 
WHERE deployment_id = 2 AND runTime != 0 AND qc_file != 2 AND 
      class_id NOT IN (58,85,88,90,91,93,94,106,109,110,111,119,121,122,124,127,128,129) 
      AND NOT (class_id = 96 AND Area > 90000) 
ORDER BY date, roi.id

ATTENTION: The subsample file in step 1 was exported from R to csv, as such, it was comma delimited while the combined file we just queried from bash is tab delimited.

The cleaning left consist in removing lines with NULL (or NA depending on your database) values because they correspond to very small cells for which the morphology could not be computed in our case (left: 110 741 171) and removing the header line. In bash, run:

grep -v "NULL" data/combined.csv > data/combined_cleaned.csv # remove lines with NULL values 
sed -i .bak '1d' data/combined_cleaned.csv # remove the header 
rm data/combined.csv */*.bak

As an example, here is how the first columns and lines of my dataset look like. As I mentioned, it has 110 741 171 images (rows) and 47 variables (columns).

To subsample, you will need the number of images per sample. The following bash code will look at the second column (filename) and count the occurence of each one of them, the second awk will just reorganize to file to have “filename,counts”. For my file the code took 45min.

time awk '{print $2}' data/combined_cleaned.csv | uniq -c | awk 'NR==1; NR > 1 {print $2","$1 | "sort"}' > data/counts.csv # the 2nd awk reorganize
sed -i .bak '1d' data/counts.csv # remove first line with the count for "filename"
head -1 data/subsp_counts_cleaned.csv | cat - data/counts.csv > data/counts_cleaned.csv # add a header based on the header of the subsampled counts (filename, count)
rm */*.bak data/counts.csv

3.2. PCA and clustering on subsampling

Subsample with the number of images you determined in 1, in my case 250.

time bash scripts/1-2_subsampling.sh -d data -p scripts -f combined_cleaned.csv \
-c counts_cleaned.csv -s date -n 250

Run the PCA and clustering on the subsampled file. Use the Rscript 3-2_pca+clustering.R.

1. Load the data (you might need to change the name of the file based on the number of images you subsampled)
2. weighted-PCA with Box-Cox transformation
3. weighted-Kmeans clustering (you might want to change the number of clusters based on what you decided)

Outputs: a csv file with the Box-Cox coefficients used for the transformation (\(\small{lambda.csv}\)), a R file with the result object of the PCA (\(\small{\textit{res_pca.rds}}\)) and a csv file with the coordinates of the centers of the morphological groups from the kmeans (\(\small{\textit{kmeans_centers.csv}}\))

3.3. Reprojection and cluster assignment

Reproject all images in the PCA space, assign a cluster number and compute the diversity indices. Use the Rscript 3-3_reprojection.R.

1. Load the data
2. Define the function to apply on each chunk of the csv file (reprojection+knn)
   • use the boxcox coefficients from the boxcox transformation of the subsample and apply them on each chunk of data
   • use the SVD matrix from the weighted-PCA run on the subsample to reproject all images in the morphological space
   • use the k-nearest-neighbour approach with the centers coordinates from the weighted-kmeans run on the subsample (k=10 neighbours) to assign a morpho group to all images
3. Run the reprojection and knn on each chunk
4. Compute the diversity indices

Outputs: a file with the coordinates in the morphological space and the cluster number for each image (\(\small{\textit{full_data_pca+kmeans.csv}}\)) and a file with the morphological diversity indices for each sample (\(\small{\textit{full_data_diversity.csv}}\))

You have your data ready for analysis!!